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Abstract 

Considering the first significant digits (noted d) in data sets of dissipation for turbulent 
flows, the probability to find a given number (d=l or 2 or... 9) would be 1/9 for an uni¬ 
form distribution. Instead the probability closely follows Newcomb-Benford’s law, namely 
P(d)=log(l-l-l/d). The discrepancies between Newcomb-Benford’s law and first-digits fre¬ 
quencies in turbulent data are analysed through Shannon’s entropy. The data sets are 
obtained with direct numerical simulations for two types of fluid flow: an isotropic case 
initialized with a Taylor-Green vortex and a channel flow. Results are in agreement with 
Newcomb-Benford’s law in nearly homogeneous cases and the discrepancies are related to 
intermittent events. Thus the scale invariance for the first significant digits, which supports 
Newcomb-Benford’s law, seems to be related to an equilibrium turbulent state, namely with 
a significant inertial range. A matlab/octave program is provided in appendix in such that 
part of the presented results can easily be replicated. 


1 Introduction 

In 1881 Newcomb [T] observed a rather strange fact: tables of logarithms in libraries tend to 
be quite dirty at the beginning and progressively cleaner throughout. This seemed indicate 
that people had more occasion to calculate with numbers beginning with 1 than with other 
digits. Newcomb concluded that the frequency of the leftmost, nonzero digit d closely follows 
the probability law: 

P(d) = logio (^1 + d = l,2,... 9 (1) 

Hence the numbers in typical statistics should have a first digit of 1 about 30% of the time, but 
a hrst digit of 9 only about 4.6% of the time (the other values can be found in tabled]). Also, 
the probability that the first significant digit is an odd number is 60%. This formula is also 
valid for the digits beyond the first, for example the distribution of the k first digits is given by 
equation dl with d = 10^“^,... 10^ — 1. 

Newcomb’s article went unnoticed until 1938 when Benford [2] has independently followed 
the same path: starting from observations about logarithm books, he deduced the same law. 
Benford also noted that this law is base-invariant, the definition m is given here with a decimal 
representation (base-10) for convenience. Until today this empirical law has been tested against 
various datasets ranging from mathematical curiosity to natural sciences [311]. 

Then there have been many attempts to rationalize this empirical law 13 El 0 El 0 do]. 
Recently Fewster m has provided an intuitive explanation based on the fact that any distribu¬ 
tion has tendency to satisfy Newcomb-Benford’s law, as long as the distribution spans several 
orders of magnitude and as long as the distribution is reasonably smooth. 
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Pietronero et al [7] have shown that Newcomb-Benford’s law is equivalent to the scale 
invariance of data set. Multiplying a data set X hy a factor A rescales the first-digits distribution 
such that: P{d[XX]) = X^P{d[X]), where d[X] denotes the first significant digit of the variable 
X and m is an exponent to be determined. The general solution of this equation takes the form 
of power-law.The probability distribution P{d), namely the sub-interval of [1, 10) occupied by 
the first digit d, is obtained by integration: 

fd+l 

P{d) = K dz 

Jd 


where K is a necessary constant to enforce the condition = 1. Newcomb-Benford’s law is 
recovered for a = 1, while a 7 ^ 1 leads to the generalized distribution found by Pietronero et 
al. [7]: 

IQi-" - 1 

The generalized Newcomb-Benford law has been applied to the distribution of leading digits 
in the prime number sequence by Tuque and Lacasa HH, the results have shown an asymp¬ 
totic evolution toward the uniform distribution (a = 0). Note that as particular case a = 1, 
Newcomb-Benford’s law satisfies a strong invariance, i.e. P(d[AX]) = P{d[X]). This result was 
first demonstrated by Hill [6]. As illustration, we consider X such that d[X] = 7 and A = 2 as 
scaling factor, thus d[AA] = 1, but the secondary digit can only be 4 or 5. With the relation 
o, it is easy to check the scale-invariance for this example: P(7) = P(14) -|- P(15). 

Self-similarity is also a well explored topic in turbulent flows. The first step was made by 
Richardson in the late 19th century with the assumption of an universal cascade process of the 
energy: the energy input at large scale is successively transferred to finer eddies. This idea 
has been refined by Kolmogorov and Obukov, the cascade is then assumed to occur in a space¬ 
filling, self-similar way. Formally there should be an unique scaling exponent for the structure 
functions S{r) such that S{Xr) = A^S'(r). Today the departures from the Kolmogorov scaling 
prediction are identified to different causes: the Reynolds number is not large enough, the flow 
is not exactly isotropic or the self-similarity assumption is not valid. The last point is related 
to the intermittent feature of the cascade process, resulting in anomalous scaling, or no unique 
scaling exponent. Further informations can be found in the book by Frisch |12j . 

In the present work, the distribution of first significant digits is used as an alternative 
statistical tool for analysing turbulent flows. Newcomb-Benford’s law is compared to data 
sets coming from numerical simulations of the Taylor-Green vortices (homogeneous case) and 
the plane Poiseuille flow (inhomogeneous case). In the next section the numerical method is 
presented and the Shannon entropy is introduced as a diagnostic tool. Results are presented in 
the following section. The validity and the possible extensions of the method are discussed in 
the last section. 



2 Method 

The momentum and mass conservation equations applied to the fluid leads to the Navier-Stokes 
equations, written in non dimensional form: 


Vu = 0, 

(3) 

clfU-I-u • Vu = —Vp-I-Pe 

the flow is assumed incompressible and the Reynolds number (Re) is the only control parameter. 
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It should be noted that the quadratic nonlinear term in the Navier-Stokes equations can be 
related to Newcomb-Benford’s law. If we assume a scale invariant data set, then the peculiar 
Newcomb-Benford’s law could be rooted in this non-linear term. Let consider the quadratic 
transformation, X —)• X^, mimicking the non linear convective term in Navier-Stokes equations 
1 The probability to find the first digit of X in the interval [d, d-l- 1) equals the probability to 
find the first digit of in [d^, (d -|- 1)^): 



dz 


leading to d = (a -|- l)/2. For any distribution, characterized by an exponent a, the quadratic 
transformation generates a distribution with an exponent d which is closer to 1 than initial a. 
As consequence the quadratic non-linear term in Navier-Stokes equations El acts as a feedback 
which force any scale invariant distribution (any a) toward Newcomb-Benford’s law (a = 1). 
However the underlying reason why scale invariance occurs is still elusive. As stated by Fewster 
m, currently any argument can explain why a universal law of nature should arise in the first 
place. 

Now we have to choose the observable which will be used to apply the statistical analysis. 
The viscous dissipation has been selected: 

^ J_ 1 ^ 

Re 2 \dxk dxi J 


In a practical point of view, the dissipation in internal or external flow fields is directly 
related to friction factor and drag coefficient respectively. In addition, for isothermal flows, the 
dissipation is the only term responsible for entropy production. Hence the dissipation can be 
considered as cornerstone for turbulent statistics, for a recent overview see Vassilicos |13] . The 
dissipation appears as a natural choice in order to exemplify the first-digit statistics in turbulent 
flows. 

In order to investigate the possible conformance of turbulent data to Newcomb-Benford’s 
law, it is necessary to quantify their discrepancy. This can be achieved through different ways. 
First the usual statistical tools used to test the goodness of fit, for example the or chi-square 
test. Secondly the search of the optimal exponent a in the generalized distribution, equation 
(©, see Luque and Lacasa m- Eventually a third choice was retained, stemming from the 
information theory. In 1948 Shannon [T4] introduced entropy as a measure of the uncertainty 
in a random variable: 

9 

H = -^P{d) logioH(d) (5) 

d=l 

This definition of entropy is associated to lack of information or to uncertainties. A message, 
in our case a chain of first digits d, occurring with probability P{d), contains the information 
/ = — log]^gP(d). Thus Shannon entropy H may be seen as the information averaged over 
the message space. That formula can be heuristically explained considering the fact that total 
information of two independent events (with respective probabilities Pi and P 2 ) is the sum 
of the two individual information, namely /(Pi, P 2 ) = /(Pi) + ^{^ 2 ), hence the logarithm. 
Moreover information is a positive quantity, hence minus sign, and the information contained 
in certain events vanishes: /(I) = 0. 

The Shannon entropy applied to the general probability distribution (equation [2]) is depicted 
in figure [Ij 

The maximum entropy is achieved when all digits are equally probable. The principle of 
maximum entropy leads to the notion of equipartition: the equal spread of probabilities over 
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Figure 1: Entropy for various distributions. Newcomb-Benford’s case, marked with a solid 
circle, corresponds to a = 1 and H = 0.8657. 

all possible states of a system. In figure [H the maximum value is H = log]^Q(9) « 0.9542. 
Conversely the minimum entropy occurs when one digit is certain, then entropy vanishes. 

The extensive property of entropy is straightforward if we consider Newcomb-Benford’s law 
for digits beyond the first. Hill [2j have shown that the probability P that d (d = 0, 1, ..., 9) is 
encountered as the n-th (n > 0) digit is: 

p„(d)= + 

The distributions of the first and second digits are not independent, however, as n increases, 
the distributions of successive digits become independent and rapidly converge to an even dis¬ 
tribution (p=0.1 for each of the ten digits). Thus, for sufficiently large n, the entropy satisfies 
the extensive property: 

H{Pni Pn+i) = H{Pn) + H{Pn+i) —)• 2 as n —oo 



3 Results 

In the following, two turbulent flows are considered. First a three-dimensional decaying and 
homogeneous turbulent flow, initialized with a Taylor-Green vortex. Then a developed wall- 
bounded turbulent flow in a plane channel. 

3.1 Homogeneous and isotropic case: the Taylor-Green vortex 

We consider a spatially periodic flow, the problem is studied by direct numerical simulation with 
256^ Fourier modes and the Reynolds number is Re = 1600. As initial condition we impose the 
single-mode Taylor-Green vortex |15[ 116] : 

u{x,y,z,t = 0) = sin(27r/3) sinx cosy cos z 

v{x,y,z,t = 0) = sin(27r/3) cosx siny cosz (6) 

w(x, y,z,t = 0) = 0 

This initial solution generates a very complex flow and evolves into turbulent flow for suf¬ 
ficiently high Reynolds number. The small scales are generated by three-dimensional vortex 
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stretching, as described in the original article by Taylor and Green m- Because of determinis¬ 
tic nature of the Navier-Stokes equations, the flow issuing from a Taylor-Green vortex presents 
microscopic structure which can be reproduced repeatedly. The results presented hereafter can 
thus be easily reproduced with the program provided in appendix. 

The time evolution of the total dissipation and the conformance of turbulent statistics with 
Benford’s law are presented in figure [2j 




Figure 2: Left: Rate of energy dissipation e as a function of time for Re = 1600. Rigth: 
Shannon entropy vs t for the Taylor-Green vortex. The dotted line stands for the Shannon 
entropy corresponding to Newcomb-Benford’s law, H = 0.8657. 

The transient increase of the dissipation, until its maximum value reached at t ~ 9, charac¬ 
terizes the nonlinear vortex stretching followed by an overall decay of the flow by viscous effect. 
For further details see Brachet et al. |16j . Shannon’s entropy slightly oscillates around the 
value expected from Newcomb-Benford’s law, namely H = 0.8657. The first-digit frequencies 
are close to Newcomb-Benford’s distribution, which is confirmed by a direct comparison shown 
in figure [3l 
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Figure 3: Gomparisons between Newcomb-Benford’s law (filled bar) and the probability distri¬ 
bution of the first significant digits of dissipation for the Taylor-Green vortex (empty bar), at 
t = 10 (left) and t = 20 (right). 
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3.2 Non-homogeneous case: the plane Poiseuille flow 

Next consider a turbulent plane Poiseuille flow, namely a channel flow driven by a stream- 
wise pressure gradient, at two different Reynolds number: Re = 2800 and Re = 6880. The 
Reynolds number is defined with the half-height of the channel and the bulk velocity. The 
no-slip boundary condition on the walls induces a different scaling based on the friction ve¬ 
locity Ur and kinematic viscosity v. Hence the wall coordinate can be expressed in wall unit: 
y+ = y vjur- The two Reynolds number values based on the skin friction velocities become 
Rer = 180 and Rer = 395, respectively. The scale separation between the inner layer and the 
outer layer increases with Reynolds number, thus it appears useful to analyse the channel flow 
for two Reynolds number values. The numerical parameters are those used in Kim et al. m 
and Moser et al. [T5] . 

Results are temporally averaged and the conformance of turbulent statistics with Benford’s 
law is represented, respectively to the wall distance expressed in outer scaling y and in wall 
units y~^, on figure [H 




Figure 4: Shannon entropy as a function of wall distance y (left) and wall distance in wall units 
y+ (right). Results for the turbulent plane Poiseuille flow, for two Reynolds number values: 
Re=2800 (thin line) and Re=6880 (thick line). The dashed line represents the entropy value 
corresponding to Newcomb-Benford’s law, namely H = 0.8657. 

Again, first-digit frequencies slightly oscillate around Newcomb-Benford’s distribution, ex¬ 
cept in the near-wall region where the oscillation amplitudes are stronger. In fact, Toschi et al. 
[Ej have observed maximum intermittency effects in that near-wall region, where it is known 
that the bursting phenomenon is the dominant dynamical feature. In addition, Toschi et al 
US! found that for positions close to the channel centreline, scaling exponents are in substantial 
agreement with observations in homogeneous and isotropic turbulence. 

Shannon’s entropy has been used as diagnostic tool to quantify the discrepancy with Newcomb- 
Benford’s law. Moreover, this information entropy can also be seen as a numerical measure which 
describes how informative a particular probability distribution is, ranging from zero (completely 
uninformative) to Hmax = log;to 9 ~ 0.9542 (completely informative). With this point of view, 
the conclusion above is consistent with results presented by Cerbus and Goldburg [20l [2T]- 
Cerbus and Goldburg observed, from experiments in a turbulent soap film, that turbulence is 
easier to predict where a cascade exists, or equivalently, for flows with a significant inertial 
range. They also found that the spatial information-entropy density is a decreasing function of 
the Reynolds number. 
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In order to complete the previous qualitative comparisons, a quantitative comparison is 
presented in the table [T] 


Table 1: Probability distribution P{d) for the first significant digits. Distribution given by 
Newcomb-Benford’s law (N-B) and distributions extracted from dissipation data for the turbu¬ 
lent plane Poiseuille flow. 


digit 

Newcomb-Benford’s 

law 

plane Poiseuille flow 
Re = 2800 

plane Poiseuille flow 
Re = 6880 

1 

0.3010 

0.2892 

0.2995 

2 

0.1761 

0.1777 

0.1724 

3 

0.1249 

0.1297 

0.1238 

4 

0.0969 

0.1011 

0.0973 

5 

0.0792 

0.0820 

0.0804 

6 

0.0669 

0.0682 

0.0684 

7 

0.0580 

0.0579 

0.0594 

8 

0.0512 

0.0501 

0.0523 

9 

0.0458 

0.0442 

0.0466 


Various tests on the numerical parameters have been realized in order to check the robustness 
of the presented results. Results remain qualitatively unchanged after: varying the spatial or 
temporal discretization, or interpolating data on different meshes, or increasing the domain size. 
Thus presented results seem to be free from possible numerical biases. 

4 Discussion 

In conclusion first-digit statistics closely follow Newcomb-Benford’s law. This is especially true 
when an equilibrium is reached: during the decay of an homogeneous turbulent flows or in 
the far from walls region of channel flows. The scale-invariant property constituting Newcomb- 
Benford’s law seems to be satisfied in conjunction with the eddies cascade in the inertial range. 

The law of first digits is not just a mathematical curiosity, the most practical use for 
Newcomb-Benford’s law is in fraud detection [22], or the possible detection of earthquakes 
[Ij. In the framework of turbulent flows, a practical application of first-digits statistics could 
be in the reduced-order models assessment, which could be investigated through their ability 
to conserve the first-digit frequencies (see Tolle et al. 123]). 

Finally the approach described in the present article can be generalized by considering the 
general form of scientific notation: 

± V lO'" (7) 

The first significant digit d{X) is the object of the present study and the two other parameters 
(namely the sign and the exponent) can also be analysed statistically. For the dissipation, 
the sign is always positive (see Equation 0]) , nonetheless for other quantities, for example the 
velocity fluctuations, the sign-change statistics provide interesting information. Indeed the time 
scale related to the fluctuating velocity zero-crossing is approximately equal to the Taylor scale 
in turbulent flows (see Kailasnath et al. [24] and references within). The last parameter in the 
formula [7| is the power m. The distributions of the power m, for the different flows considered 
here, are shown in the figure [5] 


7 








Figure 5: Probability mass function for the exponent m in dissipation data for the Taylor-Green 
vortex at t=10 (o), Plane Poiseuille at Re=2800 (A) and Re=6880 (+). 

Despite variety of the flow configurations, the curves look very similar. However in lack of 
theoretical background (as far as I know), it is not possible to go deeper in the analysis, letting 
that as an open issue. 
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A MatLab-Octave program to compute the first-digit frequen¬ 

cies 

function P=NewcombBenford(X); 

X=abs(X); 

fd = floor(X./(10.''floor(logl0(X)))); % extract the first digit 
P = histc(fd,l:9)7length(X); % probability 


B MatLab-Octave program for the Taylor-Green vortex simu¬ 
lation 


We here consider the flow of an incompressible fluid under periodic boundary conditions with 
period 27r. The pressure term can be eliminated by the incompressibility condition, then Navier- 
Stokes equations can be written as: 


9u 

'di 


u X w — V 


p + -u^ 




where a; = V x u is the vorticity. Time marching is realized with Fourth-order Runge-Kutta 
method and alias error is removed by mode truncation (Orszag 2/3 rule). 










%% simulation of Taylor-Green Vortex 


%% - 

clear all, 

Re = 1600; 

N = 2'-8; 
dt = 5e-3; 

Tend=20; 

% - 

%% Fourier pseudo-spectral method 
X = 2V/N*[0;N-1]’; kx = [0:N/2-l 0 -N/2-M:-l]’; 
for m=l;N, for n=l:N, 

IKX(:,m,n)=i*kx; IKY(m,:,n)=i*kx; IKZ(m,n,:)=i*kx; 
end, end 

K2=-{lKXr2+lKYr2+lKZr2); 

K2p=K2; K2p(l,l,l) = l; K2p(N/2+l,K2p(:,N/2+l,:) = l; K2p(:,:,N/2+l) = l; 
Z=ones(N,N,N); 

Z = 1 - (sqrt(K2) > round(N/3)-|-l); 

Z(N/2-|-l,:,:)=0; Z(:,N/2-|-l,:)=0; Z(;,;,N/2-|-l)=0; 

E = exp(-dt*K2/Re); E2 = exp(-dt/2*K2/Re); 

%% Newcomb-Benford’s law 
NBL=loglO(l-Hl./[l:9]); 

SNB(l)=-sum(NBL.*loglO(NBL)); 

%% initial condition 

[Xs,Ys, Zs] = meshgrid(x,x,x); %% 3D grid 

uf = fftn( 2/sqrt(3)*sin( 2V/3)*sin(Xs).*cos(Ys).*cos(Zs) ); 

vf = fftn( 2/sqrt(3)*sin(-2*pi/3)*cos(Xs).*sin(Ys).*cos(Zs) ); 

wf = zeros(N,N,N); 

clear Xs Ys Zs 

for k=l:round(Tend/dt) %%=========================== 

u=real(ifftn(uf)); v=real(ifftn(vf)); w=real(ifftn(wf)); 

nltu=Z.*fftn(v.*real(ifftn(IKX.*vf-IKY.*uf))-w.*real(ifftn(IKZ.*uf-IKX.*wf))); 

nltv=Z.*fftn(w.*real(ifftn(IKY.*wf-IKZ.*vf))-u.*real(ifftn(IKX.*vf-IKY.*uf))); 

nltw=Z.*fftn(u.*real(ifftn(IKZ.*uf-IKX.*wf))-v.*real(ifftn(IKY.*wf-IKZ.*vf))); 

pf=-(IKX.*nltu-MKY.*nltv-MKZ.*nltw)./K2p; pf(1,1,1)=0; 

nltu=nltu-IKX.*pf; ufa = E2.*(uf -|- dt/2*nltu); 

nltv=nltv-IKY.*pf; vfa = E2.*(vf dt/2*nltv); 

nltw=nltw-IKZ.*pf; wfa = E2.*(wf -|- dt/2*nltw); 

u=real(ifftn(ufa)); v=real(ifftn(vfa)); w=real(ifftn(wfa)); 

nltua=Z.*fftn(v.*real(ifftn(IKX.*vfa-IKY.*ufa))-w.*real(ifftn(IKZ.*ufa-IKX.*wfa))) 
nltva=Z.*fftn(w.*real(ifftn(IKY.*wfa-IKZ.*vfa))-u.*real(ifftn(IKX.*vfa-IKY.*ufa))) 
nltwa=Z.*fftn(u.*real(ifftn(IKZ.*ufa-IKX.*wfa))-v.*real(ifftn(IKY.*wfa-IKZ.*vfa))) 
pf=-(IKX.*nltua-kIKY.*nltva-kIKZ.*nltwa)./K2p; pf(1,1,1)=0; 
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nltua=nltua-IKX.*pf; ufb = E2.*(uf + dt/2*nltua); 
nltva=nltva-IKY.*pf; vfb = E2.*(vf + dt/2*nltva); 
nltwa=nltwa-IKZ.*pf; wfb = E2.*(wf + dt/2*nltwa); 

u=real(ifFtn(ufb)); v=real(ifftn(vfb)); w=real(ifFtn(wfb)); 

nltub=Z.*fFtn(v.*real(ifFtn(IKX.*vfb-IKY.*ufb))-w.*real(ifftn(IKZ.*ufb-IKX.*wfb))); 

nltvb=Z.*fftn(w.*real(ifftn(IKY.*wfb-IKZ.*vfb))-u.*real(ifftn(IKX.*vfb-IKY.*ufb))); 

nltwb=Z.*fFtn(u.*real(ifftn(IKZ.*ufb-IKX.*wfb))-v.*real(ifftn(IKY.*wfb-IKZ.*vfb))); 

pf=-(IKX.*nltub+IKY.*nltvb+IKZ.*nltwb)./K2p; pf(l,l,l)=0; 

nltub=nltub-IKX.*pf; ufb = E.*(ufa + dt*nltub); 

nltvb=nltvb-IKY.*pf; vfb = E.*(vfa + dt*nltvb); 

nltwb=nltwb-IKZ.*pf; wfb = E.*(wfa + dt*nltwb); 

u=real(ifftn(ufb)); v=real(ifftn(vfb)); w=real(ifftn(wfb)); 

nltuc=Z.*fftn(v.*real(ifftn(IKX.*vfb-IKY.*ufb))-w.*real(ifftn(IKZ.*ufb-IKX.*wfb))); 

nltvc=Z.*fftn(w.*real(ifftn(IKY.*wfb-IKZ.*vfb))-u.*real(ifftn(IKX.*vfb-IKY.*ufb))); 

nltwc=Z.*fftn(u.*real(ifftn(IKZ.*ufb-IKX.*wfb))-v.*real(ifrtn(IKY.*wfb-IKZ.*vfb))); 

pf=-(IKX.*nltuc+IKY.*nltvc+IKZ.*nltwc)./K2p; pf(l,l,l)=0; 

nltuc=nltuc-IKX. *pf; 

nltvc=nltvc-IKY.*pf; 

nltwc=nltwc-IKZ.*pf; 

uf = E.*(uf + dt/ 6 *(nltu + 2*(nltua+nltub) + nltuc)); 
vf = E.*(vf + dt/ 6 *(nltv + 2*(nltva+nltvb) + nltvc)); 
wf = E.*(wf + dt/ 6 *(nltw + 2*(nltwa+nltwb) + nltwc)); 


%% plot result 
time(k) = k*dt; 

EPS=(2*(real(ifftn(IKX.*uf)).-2+real(ifftn(IKY.*vf)).-2+ ... 
real(ifftn(IKZ.*wf)).'-2)+real(ifftn(IKY.*uf + IKX.*vf)).-2+ ... 
real(ifftn(IKZ.*vf + IKY.*wf)).'-2+real(ifftn(IKX.*wf + IKZ.*uf)).*2)/Re; 

Diss(k) = mean(mean(mean(EPS))); 

X = reshape(EPS,N''3,l); 

fd = floor(X./(10.~floor(logl0(X)))); % extract the first digit 
stat(:,k) = histc(fd,l:9)’/length(X); % compute the probability 
H(k)=-sum(stat(;,k).*loglO(stat(:,k))); 

div = max(max(max(abs( real(ifftn(IKX.*uf+IKY.*vf+IKZ.*wf)) )))); 

disp([’ time=’,num2str(time(k)),’ divergence=’,num2str(div),’ dissipation=’,num2str(Diss(k))]) 
subplot(l,2,l), plot (time,Diss), xlabel t, ylabel dissipation 
subplot(l,2,2), plot (time,H,[0 time(k)],SNB*[l l],’k’), xlabel t, ylabel H 
drawnow, 

end % TIME STEPPING ================================= 
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